This code is written for practice of data analyses of single cell sequencing data.

Data

This data is taken from 10X genomics website. This dataset consists of four dissociated lung tumours namely Lung cancer, Squamous cell carcinoma, Adenocarcinoma, and NSCLC. More information about this dataset can be found here. The code is adapted from the Seurat - Guided Clustering Tutorial which goes through the basic pre-processing of the single cell RNA sequencing data followed by clustering.

Installing the libraries

#setwd("~/Desktop/Bioinformatics/github/research-singlecellaRNAsequencing/code/")
library(Seurat)
library(tidyverse)

Loading the data set

lc.sparse.m <- Read10X_h5(filename = "../data/SC3_v3_NextGem_DI_CellPlex_CSP_DTC_Sorted_30K_Multiplex_count_raw_feature_bc_matrix.h5")
str(lc.sparse.m)
## List of 3
##  $ Gene Expression     :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
##   .. ..@ i       : int [1:91374183] 5374 8468 14933 975 11755 19623 11755 33502 31575 4311 ...
##   .. ..@ p       : int [1:3897406] 0 1 1 3 5 5 5 6 7 7 ...
##   .. ..@ Dim     : int [1:2] 36601 3897405
##   .. ..@ Dimnames:List of 2
##   .. .. ..$ : chr [1:36601] "MIR1302-2HG" "FAM138A" "OR4F5" "AL627309.1" ...
##   .. .. ..$ : chr [1:3897405] "AAACCCAAGAAACACT-1" "AAACCCAAGAAACCAT-1" "AAACCCAAGAAACCCA-1" "AAACCCAAGAAACCCG-1" ...
##   .. ..@ x       : num [1:91374183] 1 1 1 1 1 1 1 1 1 1 ...
##   .. ..@ factors : list()
##  $ Antibody Capture    :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
##   .. ..@ i       : int [1:2033572] 1 3 3 0 9 0 3 1 1 13 ...
##   .. ..@ p       : int [1:3897406] 0 0 0 0 0 0 1 1 1 1 ...
##   .. ..@ Dim     : int [1:2] 17 3897405
##   .. ..@ Dimnames:List of 2
##   .. .. ..$ : chr [1:17] "CD3" "CD4.1" "CD8a" "CD14.1" ...
##   .. .. ..$ : chr [1:3897405] "AAACCCAAGAAACACT-1" "AAACCCAAGAAACCAT-1" "AAACCCAAGAAACCCA-1" "AAACCCAAGAAACCCG-1" ...
##   .. ..@ x       : num [1:2033572] 1 1 1 1 1 1 1 1 1 1 ...
##   .. ..@ factors : list()
##  $ Multiplexing Capture:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
##   .. ..@ i       : int [1:3726686] 4 4 3 4 2 4 3 3 4 3 ...
##   .. ..@ p       : int [1:3897406] 0 1 2 2 4 4 4 5 6 6 ...
##   .. ..@ Dim     : int [1:2] 12 3897405
##   .. ..@ Dimnames:List of 2
##   .. .. ..$ : chr [1:12] "CMO301" "CMO302" "CMO303" "CMO304" ...
##   .. .. ..$ : chr [1:3897405] "AAACCCAAGAAACACT-1" "AAACCCAAGAAACCAT-1" "AAACCCAAGAAACCCA-1" "AAACCCAAGAAACCCG-1" ...
##   .. ..@ x       : num [1:3726686] 1 1 2 2 1 1 1 1 1 1 ...
##   .. ..@ factors : list()
counts <- lc.sparse.m$`Gene Expression`

Initialize Seurat Object

The data used here is raw (not normalized).

lcanSeuratObj <- CreateSeuratObject(counts = counts, project = "Lung Cancer", min.cells = 3, min.features = 200)
str(lcanSeuratObj)
## Formal class 'Seurat' [package "SeuratObject"] with 13 slots
##   ..@ assays      :List of 1
##   .. ..$ RNA:Formal class 'Assay5' [package "SeuratObject"] with 8 slots
##   .. .. .. ..@ layers    :List of 1
##   .. .. .. .. ..$ counts:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
##   .. .. .. .. .. .. ..@ i       : int [1:82542998] 72 136 172 180 391 412 694 699 765 803 ...
##   .. .. .. .. .. .. ..@ p       : int [1:65597] 0 208 1490 3084 3315 3708 3992 4229 4478 4936 ...
##   .. .. .. .. .. .. ..@ Dim     : int [1:2] 28491 65596
##   .. .. .. .. .. .. ..@ Dimnames:List of 2
##   .. .. .. .. .. .. .. ..$ : NULL
##   .. .. .. .. .. .. .. ..$ : NULL
##   .. .. .. .. .. .. ..@ x       : num [1:82542998] 1 1 1 4 1 1 1 2 1 1 ...
##   .. .. .. .. .. .. ..@ factors : list()
##   .. .. .. ..@ cells     :Formal class 'LogMap' [package "SeuratObject"] with 1 slot
##   .. .. .. .. .. ..@ .Data: logi [1:65596, 1] TRUE TRUE TRUE TRUE TRUE TRUE ...
##   .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. .. .. .. ..$ : chr [1:65596] "AAACCCAAGAGTGAAG-1" "AAACCCAAGATAGGGA-1" "AAACCCAAGATCACCT-1" "AAACCCAAGATCCCGC-1" ...
##   .. .. .. .. .. .. .. ..$ : chr "counts"
##   .. .. .. .. .. ..$ dim     : int [1:2] 65596 1
##   .. .. .. .. .. ..$ dimnames:List of 2
##   .. .. .. .. .. .. ..$ : chr [1:65596] "AAACCCAAGAGTGAAG-1" "AAACCCAAGATAGGGA-1" "AAACCCAAGATCACCT-1" "AAACCCAAGATCCCGC-1" ...
##   .. .. .. .. .. .. ..$ : chr "counts"
##   .. .. .. ..@ features  :Formal class 'LogMap' [package "SeuratObject"] with 1 slot
##   .. .. .. .. .. ..@ .Data: logi [1:28491, 1] TRUE TRUE TRUE TRUE TRUE TRUE ...
##   .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. .. .. .. ..$ : chr [1:28491] "AL627309.1" "AL627309.3" "AL627309.5" "AL627309.4" ...
##   .. .. .. .. .. .. .. ..$ : chr "counts"
##   .. .. .. .. .. ..$ dim     : int [1:2] 28491 1
##   .. .. .. .. .. ..$ dimnames:List of 2
##   .. .. .. .. .. .. ..$ : chr [1:28491] "AL627309.1" "AL627309.3" "AL627309.5" "AL627309.4" ...
##   .. .. .. .. .. .. ..$ : chr "counts"
##   .. .. .. ..@ default   : int 1
##   .. .. .. ..@ assay.orig: chr(0) 
##   .. .. .. ..@ meta.data :'data.frame':  28491 obs. of  0 variables
##   .. .. .. ..@ misc      :List of 1
##   .. .. .. .. ..$ calcN: logi TRUE
##   .. .. .. ..@ key       : chr "rna_"
##   ..@ meta.data   :'data.frame': 65596 obs. of  3 variables:
##   .. ..$ orig.ident  : Factor w/ 1 level "Lung Cancer": 1 1 1 1 1 1 1 1 1 1 ...
##   .. ..$ nCount_RNA  : num [1:65596] 252 3651 4755 285 609 ...
##   .. ..$ nFeature_RNA: int [1:65596] 208 1282 1594 231 393 284 237 249 458 250 ...
##   ..@ active.assay: chr "RNA"
##   ..@ active.ident: Factor w/ 1 level "Lung Cancer": 1 1 1 1 1 1 1 1 1 1 ...
##   .. ..- attr(*, "names")= chr [1:65596] "AAACCCAAGAGTGAAG-1" "AAACCCAAGATAGGGA-1" "AAACCCAAGATCACCT-1" "AAACCCAAGATCCCGC-1" ...
##   ..@ graphs      : list()
##   ..@ neighbors   : list()
##   ..@ reductions  : list()
##   ..@ images      : list()
##   ..@ project.name: chr "Lung Cancer"
##   ..@ misc        : list()
##   ..@ version     :Classes 'package_version', 'numeric_version'  hidden list of 1
##   .. ..$ : int [1:3] 5 0 1
##   ..@ commands    : list()
##   ..@ tools       : list()
# lcanSeuratObj1 <- CreateSeuratObject(counts = counts, project = "Lung Cancer", min.cells = 3, min.features = 200)
# str(lcanSeuratObj1)
lcanSeuratObj
## An object of class Seurat 
## 28491 features across 65596 samples within 1 assay 
## Active assay: RNA (28491 features, 0 variable features)
##  1 layer present: counts

Processing

Before we perform any statistical analyses on the data, we need to pre-process the data. The steps for doing so are as follows:

1. Quality Control

% Mitochondrial reads: In dying cells, we see a higher percentage of mitochondrial gene contamination. Therefore we use this feature to remove the reads that came from dying cell data.

lcanSeuratObj[["percent.mt"]] <- PercentageFeatureSet(lcanSeuratObj, pattern = "^MT-")
# View(lcanSeuratObj@meta.data)
VlnPlot(lcanSeuratObj, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

In the first plot we see a significant number of cells having different number of genes and is spread over a spectrum. We can see from the second plot that there are many cells from which a high number of reads were detected. In third plot we also see a high number of cells where high activity of mitochondrial genes was detected. These cells will be filtered out as they are not of good quality.

FeatureScatter(lcanSeuratObj, feature1 = "nCount_RNA", feature2 = "nFeature_RNA") +
  geom_smooth(method = "lm")

Here we see most of the cells follow the trend line. We don’t many cells that have multiple copies of same transcripts as we do not see many cells in lower right corner. We also don’t have cells were many genes are caught but not sequenced enough as we don’t see any cells in top left corner. From this graph, we can safely infer that overall, our data is of good quality. We will still filter out cells whose mitochondrial gene expression is detected above 5%.

2. Filtering

lcanSeuratObj <- subset(lcanSeuratObj, subset = nFeature_RNA > 235 & nFeature_RNA < 2500
                        & percent.mt < 5)
lcanSeuratObj
## An object of class Seurat 
## 28491 features across 25629 samples within 1 assay 
## Active assay: RNA (28491 features, 0 variable features)
##  1 layer present: counts

3. Normalizing data

Here we normalization method used is LogNormalize where Feature counts for each cell are divided by the total counts for that cell and multiplied by the scale factor. This is then natural-log transformed using log1p. The scale factor is 10000.

lcanSeuratObj <- NormalizeData(lcanSeuratObj)
# The normalization used here is standard. However parameters can be changed as below.
# lcanSeuratObj <- NormalizeData(lcanSeuratObj, normalization.method = "LogNormalize", scale.factor = 10000)

Identifying highly variable features:

Here we select the features that have high variation from cell to cell. These genes are expressed highly in some cells and lowly in others. These are our genes of interest for downstream analyses and help to point to biological signals in single cell data sets.

lcanSeuratObj <- FindVariableFeatures(lcanSeuratObj, selection.method = "vst", nfeatures = 2000)

Identifying 10 most highly variable genes: These are the top 10 genes.

top10 <- head(VariableFeatures(lcanSeuratObj), 10)

top10
##  [1] "IGHA1"  "IGKC"   "JCHAIN" "IGHG1"  "IGHM"   "TPSAB1" "IGHG3"  "TPSB2" 
##  [9] "S100A9" "IGLC2"
plot1 <- VariableFeaturePlot(lcanSeuratObj)
LabelPoints(plot = plot1, points = top10, repel = T)

Here is a visualization of the variable features where the top 10 most variable genes are highlighted.

5. Scaling

Scaling helps to remove unwanted variation from data that can be because of technical noise like batch effects or from biological sources such as different stages in cell cycle. After scaling we can have the genes clustering because of actual biological similarities and differences.

all.genes <- rownames(lcanSeuratObj)
lcanSeuratObj <- ScaleData(lcanSeuratObj, features = VariableFeatures(lcanSeuratObj))

# str(lcanSeuratObj)

6. Perform Linear Dimensionality reduction

lcanSeuratObj <- RunPCA(lcanSeuratObj, features = VariableFeatures(object = lcanSeuratObj))

# Visualize PCA results
print(lcanSeuratObj[['pca']], dims = 1:5, nfeatures = 5)
## PC_ 1 
## Positive:  CD2, IL32, TRAC, TRBC2, LTB 
## Negative:  IGKC, FTL, IGLC3, IGHA1, IGHG3 
## PC_ 2 
## Positive:  BANK1, CD79A, MS4A1, VPREB3, TNFRSF13C 
## Negative:  IL32, CD2, ANXA1, IL7R, GZMA 
## PC_ 3 
## Positive:  MCEMP1, CD68, S100A8, SMIM25, TYROBP 
## Negative:  IGKC, IGLC3, IGLC2, IGHA1, IGHG3 
## PC_ 4 
## Positive:  TNFRSF4, CTLA4, ICA1, MAGEH1, TNFRSF18 
## Negative:  NKG7, KLRD1, CCL5, CTSW, GZMH 
## PC_ 5 
## Positive:  IL7R, CD40LG, ANXA1, LINC00892, GPR171 
## Negative:  LAYN, TIGIT, TNFRSF9, LAIR2, CD27

Visualization dot plot

VizDimLoadings(lcanSeuratObj, dims = 1:2, reduction = "pca")

Dimplot

DimPlot(lcanSeuratObj, reduction = "pca")

Visualization DimHeatmap

DimHeatmap(lcanSeuratObj, dims = 1:15, cells = 500, balanced = T)

In this heatmap, we can see which features have maximum heterogeneity. Focussing on those genes in downstream process gives us a maximum sense of biological activity happening. Here we gave command to plot the 500 cells based on the PCA scores. The algorithm plots the readings with highest and lowest scores. We do see maximum of heterogeneity in PC2, PC3 and PC5; and minimum heterogeneity in PC13, PC14, and PC15.

# DoHeatmap(
#   lcanSeuratObj,
#   features = all.genes,
#   cells = 500,
#   group.by = "ident",
#   group.bar = TRUE,
#   disp.min = -2.5,
#   slot = "scale.data",
#   label = TRUE,
#   size = 5.5,
#   hjust = 0,
#   angle = 45,
#   raster = TRUE,
#   draw.lines = TRUE,
#   group.bar.height = 0.02,
#   combine = TRUE
# )

Elbow plot

ElbowPlot(lcanSeuratObj)

This is an elbow plot. Here the the elbowing is seen at 6 which means that we can get most relevant heterogeneity in the first 6 PCAs and the rest are not worth investigating further. For this parameter, it is advisable to err on higher end as not doing so adversely affects the downstream results.

7. Clustering

lcanSeuratObj <- FindNeighbors(lcanSeuratObj, dims = 1:15)
# View(lcanSeuratObj@meta.data)
lcanSeuratObj <- FindClusters(lcanSeuratObj, resolution = c(0.1, 0.3, 0.5, 0.7, 1.0))
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 25629
## Number of edges: 693785
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9674
## Number of communities: 8
## Elapsed time: 3 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 25629
## Number of edges: 693785
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9226
## Number of communities: 11
## Elapsed time: 3 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 25629
## Number of edges: 693785
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8838
## Number of communities: 14
## Elapsed time: 3 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 25629
## Number of edges: 693785
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8512
## Number of communities: 16
## Elapsed time: 3 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 25629
## Number of edges: 693785
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8211
## Number of communities: 21
## Elapsed time: 4 seconds
#View(lcanSeuratObj@meta.data)

DimPlot(lcanSeuratObj, group.by = "RNA_snn_res.0.1", label = T)

Here we see the similar cells clustering together. As the resolution increases, we get more clusters, and better defination of each cluster.

Visualiztion Dimplot with resolution 0.3

DimPlot(lcanSeuratObj, group.by = "RNA_snn_res.0.3", label = T)

Visualization Dimplot with resolution 0.5

DimPlot(lcanSeuratObj, group.by = "RNA_snn_res.0.5", label = T)

Visualization Dimplot with resolution 0.7

DimPlot(lcanSeuratObj, group.by = "RNA_snn_res.0.7", label = T)

Visualization Dimplot with resolution 1

DimPlot(lcanSeuratObj, group.by = "RNA_snn_res.1", label = T)

Setting Identity of clusters

Non Linear Dimensionality Reduction

lcanSeuratObj <- RunUMAP(lcanSeuratObj, dims = 1:15)
DimPlot(lcanSeuratObj, reduction = "umap")

# saveRDS(lcanSeuratObj , file = "../output/lungTumour-4CMO.rds")
# pheatmap(lcanSeuratObj) 
# The goal here is to cluster the cells based on cancer type

Finding Differentially Expressed Features (Cluster Biomarkers)

Finding all the markers of cluster 2:

cluster2.markers <- FindMarkers(lcanSeuratObj, ident.1 = 2, min.pct = 0.25)
head(cluster2.markers, n=5)
##           p_val avg_log2FC pct.1 pct.2 p_val_adj
## MS4A1         0   3.295984 0.905 0.039         0
## CD79A         0   3.325429 0.899 0.037         0
## BANK1         0   4.085661 0.842 0.015         0
## CD37          0   2.593567 0.993 0.279         0
## TNFRSF13C     0   3.218323 0.711 0.020         0

Finding markers of every cluster compared to all remaining cells reporting only the positive ones.

lcan.markers <- FindAllMarkers(lcanSeuratObj, only.pos = T, min.pct = 0.25, 
                               logfc.threshold = 0.25)
lcan.markers %>% 
  group_by(cluster) %>%
  slice_max(n = 2, order_by = avg_log2FC)
## # A tibble: 16 × 7
## # Groups:   cluster [8]
##       p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene     
##       <dbl>      <dbl> <dbl> <dbl>     <dbl> <fct>   <chr>    
##  1 5.91e-15       2.99 0.182 0.272  1.69e-10 0       SPINT2   
##  2 1.89e-14       2.81 0.203 0.306  5.37e-10 0       CD24     
##  3 0              3.14 0.288 0.007  0        1       CD40LG   
##  4 0              2.93 0.464 0.039  0        1       KLRB1    
##  5 0              5.06 0.604 0.005  0        2       VPREB3   
##  6 0              4.64 0.289 0.003  0        2       IGHD     
##  7 0              5.18 0.379 0.005  0        3       CD8B     
##  8 0              4.85 0.254 0.004  0        3       LINC02446
##  9 0              2.97 0.566 0.095  0        4       S100A8   
## 10 0              2.95 0.311 0.012  0        4       MCEMP1   
## 11 2.29e- 4       4.63 0.939 0.789  1   e+ 0 5       IGHA1    
## 12 0              4.47 0.537 0.007  0        5       JSRP1    
## 13 0              8.22 0.46  0.002  0        6       LINC02812
## 14 0              7.57 0.7   0.001  0        6       SCT      
## 15 0             10.1  0.306 0      0        7       CTSG     
## 16 0              9.71 0.333 0      0        7       SVOPL
cluster0.markers <- FindMarkers(lcanSeuratObj, ident.1 = 0, logfc.threshold = 0.25, test.use = "roc", only.pos = T)

VlnPlot(lcanSeuratObj, features = c("ID3", "RPL22"))

This plot can be used to visualize particular genes of interest and it shows how it is expressed in all the clusters. Based on this plot, we can examine the individual clusters where these genes are expressed most for further testing.

VlnPlot(lcanSeuratObj, features = c("ID3", "RPL22"), slot = "counts", log = T)

This is the actual raw count, before scaling and normalizing.

FeaturePlot(lcanSeuratObj,features = c("IGHA1", "IGKC", "JCHAIN", "IGHG1", "IGHM",  "TPSAB1"))

lcan.markers %>%
  group_by(cluster) %>%
  top_n(n = 10, wt = avg_log2FC) -> top10
DoHeatmap(lcanSeuratObj, features = top10$gene) + NoLegend()

# new.cluster.ids <- c("LUNG CANCER", "SQUAMOUS CELL CARCINOMA", "ADENOCARCINOMA", "NSCLC")
# names(new.cluster.ids) <- levels(lcanSeuratObj)
# lcanSeuratObj <- RenameIdents(lcanSeuratObj, new.cluster.ids)
# DimPlot(lcanSeuratObj, reduction = "umap", label = T, pt.size = 0.5) + NoLegend()
# saveRDS(lcanSeuratObj, file = "../output/lungCancerFinal.rds")
sessionInfo()
## R version 4.3.1 (2023-06-16)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS Monterey 12.2.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/New_York
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] lubridate_1.9.3    forcats_1.0.0      stringr_1.5.0      dplyr_1.1.3       
##  [5] purrr_1.0.2        readr_2.1.4        tidyr_1.3.0        tibble_3.2.1      
##  [9] ggplot2_3.4.3      tidyverse_2.0.0    Seurat_5.0.1       SeuratObject_5.0.1
## [13] sp_2.0-0          
## 
## loaded via a namespace (and not attached):
##   [1] RColorBrewer_1.1-3     rstudioapi_0.15.0      jsonlite_1.8.7        
##   [4] magrittr_2.0.3         ggbeeswarm_0.7.2       spatstat.utils_3.0-3  
##   [7] farver_2.1.1           rmarkdown_2.25         vctrs_0.6.3           
##  [10] ROCR_1.0-11            spatstat.explore_3.2-3 htmltools_0.5.6       
##  [13] sass_0.4.7             sctransform_0.4.1      parallelly_1.36.0     
##  [16] KernSmooth_2.23-21     bslib_0.5.1            htmlwidgets_1.6.2     
##  [19] ica_1.0-3              plyr_1.8.8             plotly_4.10.3         
##  [22] zoo_1.8-12             cachem_1.0.8           igraph_1.5.1          
##  [25] mime_0.12              lifecycle_1.0.3        pkgconfig_2.0.3       
##  [28] Matrix_1.6-4           R6_2.5.1               fastmap_1.1.1         
##  [31] fitdistrplus_1.1-11    future_1.33.0          shiny_1.7.5           
##  [34] digest_0.6.33          colorspace_2.1-0       patchwork_1.1.3       
##  [37] tensor_1.5             RSpectra_0.16-1        irlba_2.3.5.1         
##  [40] labeling_0.4.3         progressr_0.14.0       fansi_1.0.4           
##  [43] spatstat.sparse_3.0-2  timechange_0.2.0       mgcv_1.8-42           
##  [46] httr_1.4.7             polyclip_1.10-6        abind_1.4-5           
##  [49] compiler_4.3.1         bit64_4.0.5            withr_2.5.1           
##  [52] fastDummies_1.7.3      MASS_7.3-60            tools_4.3.1           
##  [55] vipor_0.4.7            lmtest_0.9-40          beeswarm_0.4.0        
##  [58] httpuv_1.6.11          future.apply_1.11.0    goftest_1.2-3         
##  [61] glue_1.6.2             nlme_3.1-162           promises_1.2.1        
##  [64] grid_4.3.1             Rtsne_0.16             cluster_2.1.4         
##  [67] reshape2_1.4.4         generics_0.1.3         hdf5r_1.3.8           
##  [70] gtable_0.3.4           spatstat.data_3.0-1    tzdb_0.4.0            
##  [73] data.table_1.14.8      hms_1.1.3              utf8_1.2.3            
##  [76] spatstat.geom_3.2-5    RcppAnnoy_0.0.21       ggrepel_0.9.3         
##  [79] RANN_2.6.1             pillar_1.9.0           limma_3.58.1          
##  [82] spam_2.10-0            RcppHNSW_0.5.0         later_1.3.1           
##  [85] splines_4.3.1          lattice_0.21-8         bit_4.0.5             
##  [88] survival_3.5-5         deldir_1.0-9           tidyselect_1.2.0      
##  [91] miniUI_0.1.1.1         pbapply_1.7-2          knitr_1.44            
##  [94] gridExtra_2.3          scattermore_1.2        xfun_0.40             
##  [97] statmod_1.5.0          matrixStats_1.0.0      stringi_1.7.12        
## [100] lazyeval_0.2.2         yaml_2.3.7             evaluate_0.22         
## [103] codetools_0.2-19       cli_3.6.1              uwot_0.1.16           
## [106] xtable_1.8-4           reticulate_1.32.0      munsell_0.5.0         
## [109] jquerylib_0.1.4        Rcpp_1.0.11            globals_0.16.2        
## [112] spatstat.random_3.1-6  png_0.1-8              ggrastr_1.0.2         
## [115] parallel_4.3.1         ellipsis_0.3.2         dotCall64_1.1-1       
## [118] listenv_0.9.0          viridisLite_0.4.2      scales_1.2.1          
## [121] ggridges_0.5.4         leiden_0.4.3           rlang_1.1.1           
## [124] cowplot_1.1.1